function [Y, desc] = generateFlowFeatures( flow, N, opt )
%generateFlowFeatures - Computes a flow feature.
%
% Flow format: 93 numbers
% frame_number <x_0_0> <y_0_0> <p_0_0> <x_0_1>  <y_0_1> <p_0_1> ...
%   <x_i_j> <y_i_j> <p_i_j> ... <x_4_5><y_4_5><p_4_5> x_ave y_ave
% 
% - i indexes the body quadrant of the bounding box:
%
%                0
%               1 2
%               3 4
%
% - j indexes the angle bin of a body quadrant:
%
%               2 1
%              3 . 0
%               4 5
%
% - x_ave and y_ave specifies the average flow of the box that was
%     subtracted from each angle bin.
%
  
  KTHseed = 0;
  if (opt == KTHseed)
    % Right & left movement distribution of the head
    desc = 'flow-head-right-left-distr';
    Y = [flow(1:N, idx('p', 'head', 'right')), ...
      flow(1:N, idx('p', 'head', 'left'))];
  elseif (opt == KTHseed+1)
    % Uright & Lleft movement distribution of the head
    desc = 'flow-head-uright-lleft-distr';
    Y = [flow(1:N, idx('p', 'head', 'uright')), ...
      flow(1:N, idx('p', 'head', 'lleft'))];
  elseif (opt == KTHseed+2)
    % Lright & Uleft movement distribution of the head
    desc = 'flow-head-lright-uleft-distr';
    Y = [flow(1:N, idx('p', 'head', 'lright')), ...
      flow(1:N, idx('p', 'head', 'uleft'))];
    
  elseif (opt == KTHseed+3)
    % Right & left movement distribution of the larm
    desc = 'flow-larm-right-left-distr';
    Y = [flow(1:N, idx('p', 'larm', 'right')), ...
      flow(1:N, idx('p', 'larm', 'left'))];
  elseif (opt == KTHseed+4)
    % Uright & Lleft movement distribution of the larm
    desc = 'flow-larm-uright-lleft-distr';
    Y = [flow(1:N, idx('p', 'larm', 'uright')), ...
      flow(1:N, idx('p', 'larm', 'lleft'))];
  elseif (opt == KTHseed+5)
    % Lright & Uleft movement distribution of the larm
    desc = 'flow-larm-lright-uleft-distr';
    Y = [flow(1:N, idx('p', 'larm', 'lright')), ...
      flow(1:N, idx('p', 'larm', 'uleft'))];
    
  elseif (opt == KTHseed+6)
    % Right & left movement distribution of the rarm
    desc = 'flow-rarm-right-left-distr';
    Y = [flow(1:N, idx('p', 'rarm', 'right')), ...
      flow(1:N, idx('p', 'rarm', 'left'))];
  elseif (opt == KTHseed+7)
    % Uright & Lleft movement distribution of the rarm
    desc = 'flow-rarm-uright-lleft-distr';
    Y = [flow(1:N, idx('p', 'rarm', 'uright')), ...
      flow(1:N, idx('p', 'rarm', 'lleft'))];
  elseif (opt == KTHseed+8)
    % Lright & Uleft movement distribution of the rarm
    desc = 'flow-rarm-lright-uleft-distr';
    Y = [flow(1:N, idx('p', 'rarm', 'lright')), ...
      flow(1:N, idx('p', 'rarm', 'uleft'))];
    
  elseif (opt == KTHseed+9)
    % Right & left movement distribution of the lleg
    desc = 'flow-lleg-right-left-distr';
    Y = [flow(1:N, idx('p', 'lleg', 'right')), ...
      flow(1:N, idx('p', 'lleg', 'left'))];
  elseif (opt == KTHseed+10)
    % Uright & Lleft movement distribution of the lleg
    desc = 'flow-lleg-uright-lleft-distr';
    Y = [flow(1:N, idx('p', 'lleg', 'uright')), ...
      flow(1:N, idx('p', 'lleg', 'lleft'))];
  elseif (opt == KTHseed+11)
    % Lright & Uleft movement distribution of the lleg
    desc = 'flow-lleg-lright-uleft-distr';
    Y = [flow(1:N, idx('p', 'lleg', 'lright')), ...
      flow(1:N, idx('p', 'lleg', 'uleft'))];
    
  elseif (opt == KTHseed+12)
    % Right & left movement distribution of the rleg
    desc = 'flow-rleg-right-left-distr';
    Y = [flow(1:N, idx('p', 'rleg', 'right')), ...
      flow(1:N, idx('p', 'rleg', 'left'))];
  elseif (opt == KTHseed+13)
    % Uright & Lleft movement distribution of the rleg
    desc = 'flow-rleg-uright-lleft-distr';
    Y = [flow(1:N, idx('p', 'rleg', 'uright')), ...
      flow(1:N, idx('p', 'rleg', 'lleft'))];
  elseif (opt == KTHseed+14)
    % Lright & Uleft movement distribution of the rleg
    desc = 'flow-rleg-lright-uleft-distr';
    Y = [flow(1:N, idx('p', 'rleg', 'lright')), ...
      flow(1:N, idx('p', 'rleg', 'uleft'))];
  
    
  elseif (opt == KTHseed+15)
    % xave & yave
    desc = 'flow-xave-yave-dir';
    Y = [flow(1:N, idx('xave', '', '')), ...
      flow(1:N, idx('yave', '', ''))];
    
    
  elseif (opt == KTHseed+16)
    % Right & left movement magnitude of the head
    desc = 'flow-head-right-left-mag';
    m1 = [flow(1:N, idx('x', 'head', 'right')), ...
      flow(1:N, idx('y', 'head', 'right'))];
    m2 = [flow(1:N, idx('x', 'head', 'left')), ...
      flow(1:N, idx('y', 'head', 'left'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
  elseif (opt == KTHseed+17)
    % Uright & Lleft movement magnitude of the head
    desc = 'flow-head-uright-lleft-mag';
    m1 = [flow(1:N, idx('x', 'head', 'uright')), ...
      flow(1:N, idx('y', 'head', 'uright'))];
    m2 = [flow(1:N, idx('x', 'head', 'lleft')), ...
      flow(1:N, idx('y', 'head', 'lleft'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
  elseif (opt == KTHseed+18)
    % Lright & Uleft movement magnitude of the head
    desc = 'flow-head-lright-uleft-mag';
    m1 = [flow(1:N, idx('x', 'head', 'lright')), ...
      flow(1:N, idx('y', 'head', 'lright'))];
    m2 = [flow(1:N, idx('x', 'head', 'uleft')), ...
      flow(1:N, idx('y', 'head', 'uleft'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
    
  elseif (opt == KTHseed+19)
    % Right & left movement magnitude of the larm
    desc = 'flow-larm-right-left-mag';
    m1 = [flow(1:N, idx('x', 'larm', 'right')), ...
      flow(1:N, idx('y', 'larm', 'right'))];
    m2 = [flow(1:N, idx('x', 'larm', 'left')), ...
      flow(1:N, idx('y', 'larm', 'left'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
  elseif (opt == KTHseed+20)
    % Uright & Lleft movement magnitude of the larm
    desc = 'flow-larm-uright-lleft-mag';
    m1 = [flow(1:N, idx('x', 'larm', 'uright')), ...
      flow(1:N, idx('y', 'larm', 'uright'))];
    m2 = [flow(1:N, idx('x', 'larm', 'lleft')), ...
      flow(1:N, idx('y', 'larm', 'lleft'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
  elseif (opt == KTHseed+21)
    % Lright & Uleft movement magnitude of the larm
    desc = 'flow-larm-lright-uleft-mag';
    m1 = [flow(1:N, idx('x', 'larm', 'lright')), ...
      flow(1:N, idx('y', 'larm', 'lright'))];
    m2 = [flow(1:N, idx('x', 'larm', 'uleft')), ...
      flow(1:N, idx('y', 'larm', 'uleft'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
    
  elseif (opt == KTHseed+22)
    % Right & left movement magnitude of the rarm
    desc = 'flow-rarm-right-left-mag';
    m1 = [flow(1:N, idx('x', 'rarm', 'right')), ...
      flow(1:N, idx('y', 'rarm', 'right'))];
    m2 = [flow(1:N, idx('x', 'rarm', 'left')), ...
      flow(1:N, idx('y', 'rarm', 'left'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
  elseif (opt == KTHseed+23)
    % Uright & Lleft movement magnitude of the rarm
    desc = 'flow-rarm-uright-lleft-mag';
    m1 = [flow(1:N, idx('x', 'rarm', 'uright')), ...
      flow(1:N, idx('y', 'rarm', 'uright'))];
    m2 = [flow(1:N, idx('x', 'rarm', 'lleft')), ...
      flow(1:N, idx('y', 'rarm', 'lleft'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
  elseif (opt == KTHseed+24)
    % Lright & Uleft movement magnitude of the rarm
    desc = 'flow-rarm-lright-uleft-mag';
    m1 = [flow(1:N, idx('x', 'rarm', 'lright')), ...
      flow(1:N, idx('y', 'rarm', 'lright'))];
    m2 = [flow(1:N, idx('x', 'rarm', 'uleft')), ...
      flow(1:N, idx('y', 'rarm', 'uleft'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
    
  elseif (opt == KTHseed+25)
    % Right & left movement magnitude of the lleg
    desc = 'flow-lleg-right-left-mag';
    m1 = [flow(1:N, idx('x', 'lleg', 'right')), ...
      flow(1:N, idx('y', 'lleg', 'right'))];
    m2 = [flow(1:N, idx('x', 'lleg', 'left')), ...
      flow(1:N, idx('y', 'lleg', 'left'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
  elseif (opt == KTHseed+26)
    % Uright & Lleft movement magnitude of the lleg
    desc = 'flow-lleg-uright-lleft-mag';
    m1 = [flow(1:N, idx('x', 'lleg', 'uright')), ...
      flow(1:N, idx('y', 'lleg', 'uright'))];
    m2 = [flow(1:N, idx('x', 'lleg', 'lleft')), ...
      flow(1:N, idx('y', 'lleg', 'lleft'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
  elseif (opt == KTHseed+27)
    % Lright & Uleft movement magnitude of the lleg
    desc = 'flow-lleg-lright-uleft-mag';
    m1 = [flow(1:N, idx('x', 'lleg', 'lright')), ...
      flow(1:N, idx('y', 'lleg', 'lright'))];
    m2 = [flow(1:N, idx('x', 'lleg', 'uleft')), ...
      flow(1:N, idx('y', 'lleg', 'uleft'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
    
  elseif (opt == KTHseed+28)
    % Right & left movement magnitude of the rleg
    desc = 'flow-rleg-right-left-mag';
    m1 = [flow(1:N, idx('x', 'rleg', 'right')), ...
      flow(1:N, idx('y', 'rleg', 'right'))];
    m2 = [flow(1:N, idx('x', 'rleg', 'left')), ...
      flow(1:N, idx('y', 'rleg', 'left'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
  elseif (opt == KTHseed+29)
    % Uright & Lleft movement magnitude of the rleg
    desc = 'flow-rleg-uright-lleft-mag';
    m1 = [flow(1:N, idx('x', 'rleg', 'uright')), ...
      flow(1:N, idx('y', 'rleg', 'uright'))];
    m2 = [flow(1:N, idx('x', 'rleg', 'lleft')), ...
      flow(1:N, idx('y', 'rleg', 'lleft'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
  elseif (opt == KTHseed+30)
    % Lright & Uleft movement magnitude of the rleg
    desc = 'flow-rleg-lright-uleft-mag';
    m1 = [flow(1:N, idx('x', 'rleg', 'lright')), ...
      flow(1:N, idx('y', 'rleg', 'lright'))];
    m2 = [flow(1:N, idx('x', 'rleg', 'uleft')), ...
      flow(1:N, idx('y', 'rleg', 'uleft'))];
    Y = [sqrt(sum(abs(m1).^2,2)), sqrt(sum(abs(m2).^2,2))];
    
  end
  
end

function i = idx(type, body, angle)

  if strcmpi(type, 'xave')
    i = 92;
    return
  end
  if strcmpi(type, 'yave')
    i = 93;
    return
  end
  
  if strcmpi(body, 'head')
    i = 2;
  elseif strcmpi(body, 'larm')
    i = 20;
  elseif strcmpi(body, 'rarm')
    i = 38;
  elseif strcmpi(body, 'lleg')
    i = 56;
  elseif strcmpi(body, 'rleg')
    i = 74;
  else
    error(['Unrecognized body: ' body])
  end
  
  if strcmpi(angle, 'right')
    i = i + 0;
  elseif strcmpi(angle, 'uright')
    i = i + 3;
  elseif strcmpi(angle, 'uleft')
    i = i + 6;
  elseif strcmpi(angle, 'left')
    i = i + 9;
  elseif strcmpi(angle, 'lleft')
    i = i + 12;
  elseif strcmpi(angle, 'lright')
    i = i + 15;
  else
    error(['Unrecognized angle: ' angle])
  end
  
  if strcmpi(type, 'x')
    i = i + 0;
  elseif strcmpi(type, 'y')
    i = i + 1;
  elseif strcmpi(type, 'p')
    i = i + 2;
  else
    error(['Unrecognized type: ' type])
  end
    
end

function A = angle3D(a,b)
  A = rad2deg( atan2( sum(abs( cross(a,b) ).^2,2).^(1/2) , dot(a,b,2) ) );
end

function A = angle2D(x1, y1, x2, y2)
  A = rad2deg(atan2(x1.*y2-x2.*y1,x1.*x2+y1.*y2));
end
